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ABSTRACT 

We compute model atmospheres and emergent spectra for low field {B ^ lO^^G) 
neutron stars, using new opacity and equation of state data from the OPAL 
project. These computations, incorporating improved treatments of flux transport 
and convective stability, provide spectra for hydrogen, solar abundance and iron 
atmospheres. We compare our results to high field magnetic atmospheres, available 
only for hydrogen. An application to apparently thermal flux from the low field 
millisecond pulsar PSR J0437-4715 shows that H atmospheres fit substantially better 
than Fe models. We comment on extension to high fields and the implication of these 
results for neutron star luminosities and radii. 



Subject headings: 



1. Introduction 

Thermal radiation from the surfaces of neutron stars (NS) provides important information 
about these elusive compact objects. In particular, neutron stars' thermal histories are already 
being probed by X-ray observations (Ogelmann 1995); these data place useful constraints on the 
interior physics and the equation of state of matter at super-nuclear densities (Tsuruta 1995). 
Significant limits on heating processes such as precipitation of magnetospheric particles and 
rotational energy dissipation in the crust can also be obtained. To interpret the observed X-ray 
fluxes, it is often assumed that the NS spectrum is blackbody. However in general we expect 
the emergent flux to be reprocessed by the NS atmosphere, and Romani (1987) showed that the 
departures from blackbody emissivities in the observed bands can be quite substantial. 

In Romani (1987) opacities from the Los Alamos Opacity Library (LAOL) (Huebner et al. 
1977) were used to generate model atmospheres and emergent spectra for non-magnetic neutron 
stars with a variety of surface compositions. More recent work (Shibanov et al. 1992; Potekhin 
& Pavlov, 1993) has treated the equation of state and opacity of pure H atmospheres in strong 
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(~ lO^^G) magnetic fields. This allowed Pavlov et al. (1995) to produce atmosphere models for 
high field NS; substantial departures from blackbody spectra were found, albeit not as large as 
in the non-magnetic case. Since the atomic composition of a NS surface is quite uncertain it is 
also important to treat heavy element atmospheres. Miller &: Neuhauser (1991) have calculated 
wave functions and energies of some magnetized heavy elements, and Miller (1992) produced 
polarization averaged bound-free cross sections and computed approximate model atmospheres. 
However, separate transport of the polarization modes dramatically affects the emergent spectrum 
(Pavlov et al. 1995), while at zero field bound-bound opacity can exceed bound free by large 
factors (Iglesias, Rogers &: Wilson 1992), so these heavy clement atmospheres may not adequately 
model the emergent spectrum. Thus since detailed treatments for elements heavier than H in 
strong fields are not yet available, comparison with low-field results remains useful. 

To delimit the range over which our non-magnetic atmospheres can be applied, we must 
estimate the field required to perturb the opacity significantly at frequencies of interest. In 
the presence of strong magnetic fields, the free-free opacities for one of two polarization modes 
perpendicular to the field and for radiation parallel to the field are strongly suppressed. For 
Thompson scattering, the cross section decrease occurs below the electron cyclotron energy, 
i.e. < huB = Il-6i3i2 keV, where B = 10^^i?i2G (Canute, Lodenquai & Ruderman, 1971). 
For inverse Bremsstrahlung suppression, both and kT should be less than hvB (Pavlov & 
Panov, 1976); at typical NS temperatures of ~ lO^K, these effects would affect the ROSAT 
PSPC band (0.1-2.4keV) when B12 > 0.01. Magnetic effects on the bound-free and bound-bound 
opacities depend on the full quantum states of magnetized atoms, but some generalizations 
are possible. In the case of hydrogen, detailed calculations are available (e.g. Rosner ct al. 
1984); for the field Bq = 1.2 x lO^G, at which hvB = IPYj these authors find the ground state 
binding energy Eb is increased by ~ 40%. Thus this bound-free feature too is altered when 
huB ~ Ej. The effect vanishes rapidly at lower fields, while for large fields we find their results 
give Eb ^ (2.5B / Bq)^/^ Ry. To estimate the minimum perturbing field in general, consider the 
magnetic contribution to the Hamiltonian of a hydrogenic atom (Rosner et al. 1984): 

Hr^^^/my = Hi + H2 = 2P{m±l)+l3\l (1) 

where (3 = B/{4Bq), m is the quantum number for angular momentum in the field direction, ±1 
is for the spin, and r± (in units of Bohr radii) is the atomic radius perpendicular to B. Here Hi 
is the Zeeman effect, while H2 is sometimes referred to as the quadratic Zeeman effect. For an 
energy level Eb with principal quantum number Hp and a nucleus with charge +Z, we define the 
critical field /3c to be that which gives i?mag ~ Eb- Using r± ~ np/Z'^ and Eb ~ Z'^/ripILy we 
find that Hi ~ Eb =^ (3 ~ (Zjn^pf (for m ± 1 of order unity), and H2 ~ Eb =^ /3 ~ [Z/upY. 
When Up > Z, H2 > Hi and so the critical field scales as rip'^. This scaling is borne out by 
the results of Rosner et al. (1984) for excited states of Hydrogen. However, in NS atmosphere 
conditions only modest excitation of a species is obtained before it is further ionized, so for larger 
Z the perturbation Hi dominates in almost all cases. Its dependence implies that spectral features 
will be perturbed when hus ~ Ej, in agreement with Miller &; Neuhauser (1991), who find for 



-3- 



hydrogenic atoms of atomic number Z that the lowest energy levels follow Ez{B) = Z'^Eh{B/Z^), 
where Eh is the corresponding hydrogen energy level at field B/Z"^. They also find the helium 
ground state to follow this scaling, with a screened ^cff- Line positions, which represent a 
difference between two energy states, will deviate substantially from their zero field values and 
ordering even before huB reaches either energy, but a shell feature such as the iron L-edge at 0.7 
keV in our spectra should persist to Bu ~ 0.7/11.6 = 0.06. Though the atmospheric structure will 
be affected for hvB > kT as the peak flux is altered, gross opacity features above the thermal peak 
will persevere to higher fields. Our non-magnetic atmospheres are therefore entirely suitable for 
low-field (10^ — lO^'^ G) NS such as millisecond pulsars, and indicative of some spectral features at 
higher fields. 

Two developments suggest that a re-examination of the heavy element, non-magnetic NS 
atmospheres is timely. First, improved treatments of astrophysical opacity using the OPAL code 
(Iglesias et al. 1992; and references therein) and by the OP collaboration (Seaton et al. 1994) have 
shown that the LAOL results substantially underestimate the opacity, especially for iron-group 
elements in the ~ lO^K temperature range. Use of the new opacities already seems to resolve 
several puzzles in stellar astrophysics, including discrepancies between observations of Cepheid 
pulsations and models based on LAOL data (see Iglesias et al. 1992). The improvements to the 
opacity and equation of state (EOS) are particularly important to the NS spectrum problem, as 
the photosphere in these atmospheres forms precisely at the densities and temperatures where the 
changes to the LAOL results are largest. Secondly, recent observations have resulted in a number 
of strong (ROSAT) and potential (EUVE) detections of apparently thermal emission from nearby 
neutron stars, providing limited spectral information in the EUV/soft X-ray band. In particular, 
some low field (millisecond) pulsars have now been detected. It is clear that the next generation 
of space facilities in this band will provide significant constraints on NS thermal compositions and 
emissivities. 

Thus to provide an improved baseline of models that illustrate spectral changes with surface 
composition and are directly applicable to low field neutron stars, we have constructed neutron star 
atmospheres and emergent spectra using the new opacities and consistent EOS, for pure iron, pure 
hydrogen, and the solar abundances of Grevesse &i Noels (1993). Improvements in the radiative 
transfer are also introduced and the resulting spectra should superceed those of Romani (1987). 
Initial results of these computations have been reported in Romani, Rajagopal, Rogers & Iglesias 
(1995). In Section II we describe the calculation of the model atmosphere structures; in Section 
III we discuss convectional stability of the atmospheres, and present emergent spectra. Spectra 
for various compositions are compared to black bodies and the high field hydrogen atmospheres 
of Pavlov et al. (1995). In Section IV we provide an initial application, fitting the spectra to 
ROSAT observations of the millisecond pulsar J0437-4715. Section V describes implications of 
these results and prospects for future work. 
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2. Model Atmospheres 



2.1. 



OPAL Data 



Because the EOS treatment allows extension to high density, OPAL EOS and opacity tables 
(kindly supplied by F. Rogers and C. Iglcsias) are our primary data source for these models. 
OPAL opacity computations are based on the method of detailed configuration accounting using 
LS coupling (Iglesias, Rogers, &: Wilson, 1987), and now use full intermediate coupling for iron 
to incorporate spin-orbit interactions (Iglesias et al. 1992). Abundances of all possible ions are 
obtained from an activity expansion of the grand canonical ensemble (Rogers 1986), and tables 
of the equation of state (EOS) completely consistent with that used for the opacity calculations 
are available (Rogers, Swenson, & Iglesias, 1995). For the present models, the EOS tables provide 
pressure and other thermodynamic variables as a function of temperature and density, interpolated 
using routines supplied with the data. 

As an example, for pure iron the opacity grid was obtained for each decade of i? = p/Tg , 
where Tq = T/(10® K), from -5 to +1; with Log(T) from 4.5 to 7.5 every 0.25, and with 10^ 
photon energies linearly spaced over 0.002 < u = E^/{kBT) < 20.0. As needed, we interpolate 
logK linearly against logi? and logT and linearly against u (at constant p and Ej). To allow 
extension to high photon energies, we use additional tables for the same R and T spaced linearly 
in E^ from leV to lOkeV at leV intervals. These tables, used only when E^ > 20A;bT, allow us 
to solve the transfer equation for Ey up to 10 keV, through the lower temperatures present in the 
outer atmospheres. The data for the pure H and solar (Grevesse &; Noels 1993) abundances are 
similar, densely covering the thermal peak with extension up to 1 keV. 



We define ab initio a vector of 120 mean optical depths tt, logarithmically spaced between 
10~^ and 10^-^. These correspond to physical depths through dz = dTx/pKr, where kt is the 
total mean opacity in cm^g~^ including electron conductivity (more below). The temperature 
and density at these grid points are adjusted until both hydrostatic equilibrium and energy 
steady state (total astrophysical flux F = ctT^j/tt) are achieved at all depths. For each of the 
three compositions, we generate six model atmospheres with T^g from 10^'^^ to 10^'^, four per 
decade. For the light element mixtures, the EOS tables stop before tt = 10^'^ in the two coolest 
atmospheres, and for iron in the one coolest and two hottest. In these cases, the depth grid 
extends to tt = 10^^^; this affects only the highest E^ slightly. 

At each frequency v, the outward flux F,^ at the frequency-specific optical depth t^, is obtained 
from the Milne integral 



2.2. Model Atmosphere Construction 




= 2 / 5,(r^)i?2(r^-r,)dr^ 




(2) 
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where Sy is the radiation source function, and E2 is the second exponential integral, which has 
maximum of one at zero, with decay length of order unity (Mihalas 1978). Similarly, the intensity 
Jy is given by 

Ju{ru) = 7T / S,{tI)Ei \tI - T,\ dr^ (3) 
^ Jo 

where Ei is the first exponential integral. We assume Local Thermodynamic Equilibrium, taking 
the source of radiation at each depth to be a Planck function at the local temperature By[T). 
The atmosphere models begin with the grey opacity solution r^(r) = (3/4)Tg^g(r + <?(t)), using 
tabulated values for ^(t). In general, in a non-grey purely radiative atmosphere, at depths so large 
that radiative flux is given at all frequencies by the diffusion approximation 

4/1 dBy\ fdT\ ^ f ^ ^ 16a .dT 

the grey solution is obtained on the optical depth scale defined by the Rosseland mean opacity. In 
the dense atmosphere of a neutron star, however, energy transport by electron conduction can be 
significant. We thus define an equivalent conduction opacity Kc by analogy with equation (|^ so 
that —Ac dT/dz = Fcond = (16(Tr^)/(3/?Kc) dT/dz, where Ac is the thermal conductivity. Here 
we use Kc = (2.5 x W^/Ne) (Z'^/A) (Tj^"^ / p) cm^/g where N^, is the number of free electrons 
per AMU, A is the average atomic mass, Z = AN^ the average ionic charge, T7 is in 10^ K, and 
p is in g/cm^. Adding harmonically with the radiative opacity provides a total mean opacity 
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R 



+ ^ on whose depth scale the initial grey model provides the correct atmosphere 



solution at large r. 

We start with the grey atmosphere temperature structure, and impose hydrostatic equilibrium 
starting at the surface (P = 0). For each tt layer, using an initial density po from the layer above, 
we calculate the pressure at the base of that layer from dP/drji = g/ kr{pq,T). We then obtain 
a corrected pi(To,Po) from the EOS tables, and iterate to convergence. At each iteration, we 
extract the opacity on an 800 bin logarithmic frequency grid extending to lOkeV, saving harmonic 
means of 50 OPAL opacities sampled within the bin. This grid of values is used to produce the 
Rosseland mean opacity kr. Once convergence is reached at each depth, the final k,^ and are 
saved, and similar arithmetic bin averages used to produce the Planck mean Kp. Once hydrostatic 
equilibrium is achieved for the whole atmosphere, we solve the transfer equation (|2|) for each of 
the same 800 frequencies, first computing t^, corresponding to each of the 120 tt by integrating 
dTy = {ku/ Kfi)dTji down from the surface. At very large depth, where Sy{Ty) varies slowly, we are 
in the diffusion limit and Fy is obtained from equation the intensity J^, always comes directly 
from equation (|3|). Total flux F and intensity J are obtained by summing over the frequency grid, 
while the flux and absorption mean opacities Kp and nj are obtained by appropriately weighted 
sums over the k^,. 

The temperature structure is then corrected using the Lucy-Unsold procedure to drive 
towards the desired steady state flux: 



= 16^ 



J r '^I^ AF{T')dT' + 2AF(0) , 

up V Jo I^r[t') J Kp dTR 



(5) 
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(Mihalas, 1978) where r = tt as defined above, and AF = aT^^/i^ — Frad — -^cond- The corrected 
temperature run is then smoothed and any temperature inversions in the outer atmosphere are 
removed. The procedure is iterated to convergence {irAF/aT^Q < 0.005 at all depths). 



3. Results 




-4 -2 

Log(p[g/cm^]) 

Fig. 1. — Iron atmosphere structures. Open circle regions of the curves are convectively unstable, 
and size of circles in the field indicates the adiabatic index of the gas. 

Wc show the converged temperature vs. density runs for iron in Figure 1. Open circles 
indicate regions of convective instability, by definition where (dlogr/dlogP)atm = Aatm > = 
(dlogT/dlogP)adiabatic- The valuc of Aad for iron from the EOS is indicated by the size of the 
background circles; instability in the radiative/conductive solutions occurs chiefly where it is low, 
in ionization zones. To gauge the importance of this instability we compute the energy gain ratio 
of a rising, but radiating, bubble 



Kp^Cp 



(6) 



{e.g. Bohm-Vitense, 1992) where v is the mean bubble speed and we take the bubble mean free 
path Zmfp to be the atmospheric scale height Pdz/dP. For small F, departures from the radiative 
temperature gradient scale as 9r^/4. If the bubbles rise freely and quasi-adiabatically, then equating 
the potential energy gain to the kinetic energy gives v ~ [(7(dT/dz|rad — dr/dz|ad)/r]^^^^mfp as 
an upper limit to the bubble velocity. In this worst case picture, F is as large as ~ 5 at large 
depth in the LogT = 5.75 atmosphere. Accordingly the emergent flux at several keV, where Tj, 
is modest at large depth, will be slightly affected. However, even the lowest known magnetic 
fields for neutron stars are ~ lO^G. If B'^/Stt > pg{dT/dz\rad - dT/dzlad) ^mfp/^j which is true 
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for B > 10^ G throughout all the atmospheres we have computed, then the magnetic field will 
suppress the convection. Mass motions can then only occur on the diffusion time scale, which for 
a partly ionized plasma allows fdiff ~ 1-1 x 10*^Tg ^''^(Icm/Zmfp) cms~i^. Under these conditions 
r is no larger than 0.04 except in the coldest iron atmosphere, in which convection might cause 
departure from the radiative temperature gradient at the 10% level. Accordingly, even when the 
magnetic field is too small to affect the EOS and opacities it should strongly suppress convection 
in most cases. 
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Fig. 2. — Emergent spectra for atmospheres containing heavy elements. See text for comments. 

We calculate emergent spectra by solving the transfer equation at the edge of the atmosphere. 
To resolve lines from intermediate depths, we compute emergent spectra on a 2 x 10^ point energy 
grid, logarithmically spaced between leV and lOkeV. The spectra, binned down to 10"^ energies, 
are shown in Figure 2. Much coarser binning still is appropriate to the resolution of current NS 
observations, and as argued above, broad-band features will be stable even to moderate NS fields. 
The iron spectra all show features at the K, L and M edges, which are also present in the solar 
abundance spectra. Narrow absorption lines are formed at small physical depth, while strong 
pressure broadening affects lines from the deep layers. All three compositions are compared in 
Figure 3 (left); the hydrogen spectra show the expected hardening due to the dependence of 
the opacity, which allows flux to shift above the Wien peak. 

The two major spectral effects of magnetic field on hydrogen are visible in Figure 3 (right). 
First, the high energy hardening is not nearly as marked, since the diminished magnetic opacities 
which dominate transfer fall much less steeply than (Shibanov et al. 1992, Fig. 1). Secondly, 
the higher photoionization threshold increases absorption between ~ 0.1 and 0.3 keV for the fields 
shown; the sharp onset is smoothed out by pressure effects on the atomic initial state (Pavlov et 
al. 1995). At the lower field, this effect causes the prominent dip in the logT = 5.5 spectrum; at 
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logT = 6.0 the dip is minimal, due to the lesser abundance of atomic hydrogen. At the higher 
field, only the onset is seen in the log T = 5.5 spectrum, and it is clear this edge could dominate a 
fit. 



I ' ' ' ' I ' ' ' ' I ' ' ' ' I ' ' ' ' I 

Composition effects 
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Fig. 3. — Left - comparison of different atmospheric compositions at two T^ff- Right - comparison 
of hydrogen spectra from the present work with the magnetic atmospheres of Pavlov et al. (1995). 



4. Application 

To apply these spectral results to soft X-ray data, we convert our spectra for each composition 
into a tabulated model for the XSPEC package. Our model supplies unredshifted flux at the NS 
surface. These spectra are subjected to a specified gravitational redshift, and overall normalization 
is fit as a free parameter. Model surface fiux Fns is related to model observed flux Fobs by: 

i^obs(^7) - — jYVz) — \sw^) \TVz) ' 

where area ^em of the neutron star at distance D emits into solid angle $7; Z is the gravitational 
redshift, and the last factor due to time dilation. For full surface emission Q = 47r; if the star has 
a small emission region such as a heated polar cap wc take Vt = 2tt, noting that limb-darkening 
would tend to decrease this value while gravitational light bending near the NS surface would 
cause it to increase. 

As an example, we fit our models and magnetic hydrogen atmosphere emergent spectra for 
1.2 X 10^^ G and 4.7 x 10^^ G, (Pavlov et al, 1995), to a hypothesized thermal component in the 
X-ray emission of the nearby millisecond (P = 5.7 ms) pulsar J0437-4715 measured by ROSAT. 
Due to its low inferred surface dipole field of 8.3 x lO^G this is an ideal candidate, though its 
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Fig. 4. — Best fits of model thermal components plus off' pulse power law, to J0437-4715 ROSAT 
PSPC data in a crucial energy range. The model fits to the high-count soft peak below 0.3 keV are 
indistinguishable. Pi channels of the PSPC detector are here binned in groups of five, for display 
purposes only. 

large spindown age ~ 7.6 x 10^ yr (Johnston et al. 1993) suggests that any thermal emission 
should come from reheating. Wc have re-analyzed the September 20-21 1992 PSPC data studied 
by Becker &: Trumper (1993). This data set is in three segments, each of about 2000 seconds' 
duration, the last two about one hour apart and separated from the first by about one day. We 
use 1300 photons from a circle 2 arcmin in radius centered at the X-ray position. Arrival times at 
the spacecraft clock are barycentred using standard IRAF PROS routines, and corrected for the 
pulsar's orbital position using an updated radio ephemeris (Johnston et al. 1993; Bell, 1995). The 
nominal phase of the X-ray pulse in the last two segments aligns with the radio peak to within 
the ROSAT clock accuracy. An apparent ROSAT clock error after the first data segment required 
us to phase it separately, and make the necessary shift to align it with the other two. Based on 
the resulting light curve, we divide the data set into on- and off-pulse phases of equal duration; 
the on-pulse interval contains ~ 60% of the counts. While the clock error decreases confidence 
in absolute phasing somewhat, the measured alignment of the X-ray pulse with the radio peak 
(inferred from polarization data to be coincident with the neutron star magnetic axis; Manchester 
&; Johnston 1995) supports the picture that this pulse represents thermal emission from a 
reheated polar cap. Hypothesizing that cap emission is superimposed on an unpulsed background 
magnetospheric or plerionic in origin, we fit a power law to the off-pulse data, finding a power law 
index of 2.5, with absorption column density of 1.1 x 10^*^ cm~^, in agreement with the results 
of Becker & Triimpcr (1993). We then added each thermal candidate model to the frozen power 
law with absorption held fixed, finding the best-fit parameters to the on-pulse data for blackbody 
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spectra, atmosphere spectra for our three compositions, and spectra from the magnetic hydrogen 
atmospheres. Ah models are processed through the current PSPC response matrix for comparison 
to the pi-channel data, and given Z = 0.306 redshift for a standard 1.4M0 neutron star 10km in 
radius. Fit parameters are determined via the maximum likelihood method using the C-statistic 
on unbinned data; the for these fits are then determined from data binned to at least 20 counts 
per bin, with bin-width at least 1/3 the half-max width of the ROSAT energy response (Table 1). 
There are 12 DOF for the magnetic models, 11 for the others. Figure 4 shows the best fit models 
for black body and for each composition, smoothed through the detector response. Given the 
fitted normalization of the thermal component, we use equation (|^, the radio dispersion pulsar 
distance of 140pc, Z = 0.306 as above, and = 27r to infer the emission areas (Table 1). 



Table 1: Fits to PSR J0437-4715 PSPC data 



Model 


Teff(106K) 


Aefr(km2) 


xVdof 




BB 


1.80 


0.31 


1.43 


.205 


H 


0.60 


22 


0.95 


.605 


Fe 


1.82 


0.34 


2.5 


.002 


Sol 


1.4 


0.73 


1.5 


.166 


Mag(1.2) 


1.0 


4.3 


0.91 


.642 


Mag(4.7) 


1.0 


4.2 


0.93 


.623 



The fit for both atmospheres containing iron is significantly poorer because the strong L-edge 
(see Fig. 2) provides excess flux near 0.7 keV where the data show a deficit. Further, the model 
counts fall sharply above the edge where the data show an excess. The problem is exacerbated by 
redshifting the model edge to 0.6 or 0.5 keV. The BB and H fits are acceptable. Following our 
hypothesis that the thermal emission comes from re-heated polar caps, we can check the fitted 
effective areas against the cap area expected in the aligned dipole rotator model: the polar cap 
radius is rpc = r sin~^(r/ric)^''^, where r is the pulsar radius and Hc = {Pc/2tt) the light cylinder 
radius. In the case of J0437-4715, for r = 10 km we have rpc = 1.9 km, giving A ^ vrrp^. = 12 km^. 
Only the H atmosphere fits give areas close to this standard polar cap; BB and Fe fits imply 
areas ~ 40 times too small. While the magnetic H fits and areas are acceptable, the 4 x lO^G 
light cylinder field inferred from spindown torques would not allow surface B ~ 10^^ G unless the 
magnetic structure were of octupole order or higher. The resulting cap structure and thermal 
pulse would differ greatly from the dipole values inferred for PSR J0437-4715. 

5. Conclusion 

The behavior of NS spectra at energies high in the ROSAT band is crucial for temperature 
determination, as interstellar absorption often leaves only the Wien-like tail detectable. In that 



- 11 - 



case, the inferred temperature depends strongly on the composition model adopted: in most cases 
the atmospheric temperatures are substantially lower than those inferred for blackbodies, although 
high temperature Fe atmospheres may require a Tcfj slightly higher than black body. Even with 
present limited spectral information and moderate absorption this composition dependence is 
significant. As an illustration, we have drawn sample ROSAT PSPC data sets as seen with 
absorption of lO^'^ cm~^ from our redshifted {Z = 0.306) models, and fit them with black body 
spectra. Figure 5 shows the ratios of fitted to actual temperature. For the B12 = 4.7 hydrogen 
atmosphere at it's lowest two temperatures, the black body temperature is driven down by the 
photoabsorption edge (above, and Figure 3). Standard cooling curves for neutron stars (e.g. 
modified URCA processes) run about four times hotter than curves for more exotic processes (e.g. 
direct URCA, pion condensates) before photon cooling takes over (Ogelman, 1995; Tsuruta, 1995) 
at t ~ lO^y. Black body fits to thermal spectra from Ogelman's four initial cooling candidates 
place them quite near the standard curve, but Figure 5 demonstrates that even magnetic light 
element atmospheres could move their correct positions close enough to the cooler exotic process 
cooling to confuse matters considerably. 



Black Body Fits to Model Data 




5.5 6 6.5 

Log(Te„) 

Fig. 5. — Fractional temperature error incurred when simulated 3000-count ROSAT PSPC 
observations of our model spectra are fitted with a black body. 

Our emergent spectra are available in electronic form for the reader interested in determining 
flux implications for a particular band or detector, cither as text files or as XSPEC models. For the 
pulsar J0437-4715, enough spectral information exists to break the composition/Tcfj degeneracy, 
favoring a pure hydrogen atmosphere quite strongly. The implied polar cap area is in much better 
agreement with theoretical expectation than that inferred from a blackbody fit. Future missions 
(e.g. AXAF &; XMM) should be able to make similar measurements of other millisecond pulsars. 
Interestingly, non-magnetic fits to the 'thermal' emission of young, high-field pulsars seem in many 
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cases to favor H atmosphere models to blackbody spectra (Ogelman 1995 and references therein). 
While this agrees with the the compositional inference from PSR J0437-4715, comparison with 
heavy element magnetic atmospheres is needed to substantiate this result. A correct treatment of 
the detailed spectrum of such high B neutron star atmospheres will require extensive improvements 
to present atomic structure and opacity computations, although approximate models following 
broad-band features may now be feasible (work in progress). For younger pulsars, such models 
together with data from the next generation of X-ray satellites should enable a serious study of 
NS surface conditions, including composition, magnetic field and even redshift, with important 
implications for NS evolution and the EOS of matter at very high densities. 

We are indebted to Forrest Rogers and Carlos Iglesias for computation of the OPAL data 

for the Fe models and assistance with their interpretation. Wc also thank G. G. Pavlov for 
providing high-field hydrogen spectra, I.-A. Yadigaroglu for assistance with the J0437 data, and 
the referee for a careful reading. RWR was supported in part by an Alfred P. Sloan fellowship 
and NASA grant NAGW-2963, and MR in part by a fellowship from the National Science and 
Engineering Research Council of Canada. This research has made use of data obtained through 
the High Energy Astrophysics Science Archive Research Center Online Service, provided by the 
NASA-Goddard Space Flight Center. 
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